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I. OUTLINE 

A set of k mutually unbiased bases (MUB) in is a set of orthonormal bases whose basis 
vectors obey the following relations: 



where a, [3 = 1 , . . . , k and i, j = 1, . . . , N. 

We will describe a few approaches to the notorious problem of proving the (non)-existence of 
four mutually unbiased bases in dimension 6. These will include the notions of Grassmannian 
distance, quadratic matrix programming, semidefinite relaxations to polynomial programming, as 
well as various tools from algebraic geometry. 

H. GRASSMANNIAN DISTANCE 

A unit ket \e) in C 6 is identified with the density operator \e) (e\ in the real affine space of unit- 
trace hermitian operators acting on C 6 . Choosing the completely mixed state as the origin, this 
becomes a vector space of real dimension 35, which can also be viewed as the space of traceless 
hermitian operators on C 6 . Thus \e) corresponds to the real vector e = \e) (e\ - gl. The inner 
product between two vectors is defined via the traceless matrices that they are paired with, that is, 
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if m, <-> Mj - \ I, ?' = 1,2, then mi m 2 = ^tr {(Mi - \l)(M 2 - |l)}. Therefore, we can view the set 
of unit ket vectors as an embedded subset of K. 35 . 

An orthonormal basis of ket vectors ik,)},=i,...,6 corresponds to the vectors {e,} ;= i ... g, which 
together span a 5-dimensional subspace tlQ. Furthermore, the condition of mutual unbiasedness, 
\{ e i \fj)\ = e>' becomes an orthogonality condition, e,- • f ; = 0, between the subspaces representing 
the bases Wei)}, {\fj)}. The Grassmannian of 5-planes in K 35 can be made into a metric space as 
follows. Let IL denote the orthogonal projector onto a 5-plane. Then the function D(ni,n 2 ) 2 = 
^tr {(IIi - II2) 2 } is the desired distance function. Note that D 2 e [0, 5], with the maximal distance 
attained iff 111 and U 2 are mutually orthogonal. We can extend this notion of distance to an average 
distance, defined as a function of four rank-5 projectors: D(IIi, n 2 , n 3 , IT4) 2 = ^ Yui<j D(Hu II;) 2 . 
We also have D 2 e [0, 5], with the maximal average distance attained iff each pair of projectors is 
mutually orthogonal. 

Since an orthonormal basis in C 6 can be represented by some rank-5 projector on K 35 , we 
can study the existence of mutually unbiased bases by looking those projectors. The idea is to 
maximize the function D 2 over quartets of rank-5 projectors representing bases in C 6 . A global 
maxima that is strictly smaller than 5 then suffices to prove the non-existence of four mutually 
unbiased bases in C 6 . There are, however, major problems with this approach. 

Firstly, there is the troublesome constraint that the IT, are rank-5 projectors which come from 
orthonormal bases in C 6 . Presumably, one parameterizes the II, by regarding them as elements of 
the vector space of symmetric 35 x 35 matrices (or simply as matrices with an additional sym- 
metry constraint). Then, the rank-5 projector property is imposed by the constraints IT 2 - IL = 
and tr{IL} = 5. As a first simplification, we ignore the requirement that the IT, correspond to 
orthonormal bases. The objective function being maximized is quadratic in the parameters speci- 
fying the matrices II, . Altogether, we have at least a quadratically-constrained quadratic program 
(QCQP), which in the general case is A^P-hard. Occasionally, if the quadratic forms involved are 
definite, one can use Schur complements to turn the QCQP into a semidefinite program, whose 
global maxima can of course be found. This is, unfortunately, not the case here. We might ask 
for something less, such as an upper bound for the global maxima, rather than its actual value. 
However, the bounds obtained so far have been trivial. One can also look at the Lagrange dual 
of this (primal) optimization problem, the reason being that the dual gives upper bounds to the 
primal, and furthermore, the dual is always convex. Strangely enough, the dual problem gives 
minimally trivial upper bounds, which roughly speaking, indicates that the non-convexity of the 
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primal is crucial and must be taken into account if we wish to draw non-trivial conclusions. And 
we have not even checked that the projectors refer to bases. We cannot even hope to relax this last 
constraint, because we would be left with the equivalent geometrical problem of packing 5-planes 
orthogonally in R 35 , for which it is trivially true that a maximum of exactly seven such planes can 
be fitted. 

Of course, one can start directly from the bases of ket vectors and then build the corre- 
sponding rank-5 projectors. However, the average-distance function D 2 becomes quartic, which 
as an optimization problem is even more difficult than a quadratic one. 



HI. QUADRATIC MATRIX PROGRAMMING AND DATTORRO'S CONVEX ITERATION 

Some consolation can be derived from the fact that a number of techniques exist for handling 
QCQP. In particular, semidefinite relaxations of QCQP have been studied for some time. Of 
interest here is a variant of this idea, which is called quadratic matrix programming Q2|] (QMP). 
This refers to nonconvex quadratic optimization problems of the following form: 



minimize tr (x r A xj + 2tr \bIx\ + c 

subject to tr jx T A ; x} + 2tr {bJx} + c t = 0, / = 1, . . . , k, 

where A t e R' IX ", B t e R nxr , c t e R, i = 0, . . . , k. Note that the objective and constraint functions 
are quadratic matrix functions of order r. We shall use the abbreviation fi(X) = tr jx r A,xJ + 
2tr {-Bf -X] + Ci, and call fi homogeneous if B t = nxr , Q = 0. 

The construction of a semidefinite relaxation to the QMP © begins by the process of homog- 
enization. For /; as given above, the corresponding homogenized quadratic matrix function ff 1 is 
defined by 



f*(Y; Z) = tr jy r A,y} + 2tr [z t bJy] + -f tr [z T z] , Y e R nxr , Z e R rxr , (3) 
which is homogeneous, and can be represented by the matrix 

flY;Z) = tr {(Y,ZfM(fi)(Y,Z)}. (4) 
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The homogenized version of © reads 

minimize f"(Y;Z) 

subject to /f (T; Z) = 0,i = l,...,k, (5) 
^ ij (Y;Z) = 26 ij ,l <i<j<r, 

where i{/,j(Y;Z) = tr \^Z T (E\. + E^z} with £T. defined as the r x r matrix with zeroes everywhere 
except for the (i, j)-th entry, which is 1 . The importance of © lies in the fact that it is solvable 
precisely when the original QMP © is solvable, and in that case, the optimal values of both 
problems are equal Ofl. 

With the homogenized problem at hand, we can form a semidefinite relaxation as follows. We 
write (Y;Z) as the single matrix variable W e R.(' !+ '") x '" 5 so that © turns into 



minimize tr {M(f )WW 7 
subject to tr [M(fi)WW T ) = 0, / = 1, . . . , k, (6) 
tr{Ar 7 WW r } = 26ij, \<i<j<r, 



where 

(7) 



N,j = 



F' + F 



and the cyclic property of the trace has been been used. Now, observe that U = WW T is a 
symmetric, positive semidefinite (n + r) x (n + r) matrix, so we have, equivalently, the optimization 
problem 

minimize tr{M(/ )C/} 

subject to tr {M(fj)U} = 0, / = 1, . . . , k, 

tv{N ij u} = 26 ij ,l<i<j<r, (8) 



U>0 

rank(i7) < r. 

Notice that, apart from the rank constraint, © is a semidefinite program. 

Let us try to relate QMP to the MUB existence problem. We have already seen, from the 
Grassmannian distance approach, that the basis vectors making up the candidate MUBs should be 
specified individually, because it is not clear how to ensure that a rank-5 projector in R 35 corre- 
sponds to a basis in C 6 . If we choose to specify ket vectors, the MU conditions invariably become 
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quartic. Therefore, we choose a compromise, which is to specify rank-1 density matrices. Note 
that two such density matrices p\ = \e) (e\ ,p 2 = |/) (f\ are mutually unbiased iff their Hilbert- 
Schmidt inner product tr {pip 2 } = (e \ f) (f\e) = g. 

In order to reduce the number of variables in the QMP that we are about to build, we can 
consider MU constellations yfl rather than full bases. In the density operator picture, these are sets 



1 i if a + p. 



of pure states p? which obey 



\a=\,...,A 
l/=l,...,6 

these pure states also obey the MU conditions ©. A MU constellation is labelled by a set of 



[ \a=l,...,4 , 

A necessary condition for a set of four MUBs [pf 6 to exist in C is that every subset of 



numbers {a\,a 2 , ■ . . ,0.^}^, which indicates that there are k groups of pure states, and that the i- 
th group comprises a, pairwise orthogonal pure states. One of the smallest constellation that is 
not known to exist is {5,3,3,3)6. By applying a global unitary transformation (or a change in 
basis), we may assume that the five pure states in the first group are given by the diagonal matrices 
E%, i = 1, . . . , 5. This leaves us to specify three groups of three pure states jpf }. , which we 
accomplish by specifying their real and imaginary parts (so eighteen 6x6 real matrices have to 
be specified). For instance, p" = C" + iD". The matrix X appearing in © will then be the vertical 
concatenation of these 18 matrices. 

It remains to verify that the objective function as well as all the constraints that we need to 
impose are in fact quadratic matrix functions. We will try to set up a feasibility problem, so the 
objective function is just the zero function. The first constraint will be that C" + iD" is hermitian, 
which is equivalent to saying that Cf is symmetric and is antisymmetric. Since the symmetric 
and antisymmetric subspaces of R 6x6 are mutually orthogonal, we just require that the components 
of Cf (resp. D") in a basis of antisymmetric (resp. symmetric) matrices vanish. These conditions 
are of the form tr {C"5 ant i symm } = and tr ^D" 5 symm } = 0, which are quadratic (linear, in fact). 

Next, we ensure that p" = C" + iD" is a rank-1 projector by imposing tr {p^} -1=0 and 
ip") 2 - p" = 0. It is straightforward to see that these conditions are quadratic in Cf and D" . 
Finally, the inner product between p a and pP. reads 

tr {pfp^} = tr {(C? + iDf)(dj + ilfj)} = tr jcfC^ - D^} , (10) 



which is again quadratic. Therefore, the MU conditions can be imposed via quadratic matrix 
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functions. 

Following the prescription described above, one arrives at a SDP with a rank constraint 
(rank(£/) < 6), which is equivalent to the original feasibility problem. In order to obtain a certifi- 
cate of infeasibilty, we can do the following. Convert one of the orthogonality conditions in the 
MU constraints, say tr jpjpi,} = 0, to an objective function f . We see that fo is non-negative, and 
attains zero precisely when orthogonality, along with all the other constraints are fulfilled. That 
is, min(/ ) = iff the {5,3,3,3)6 constellation exists. Therefore, if we can show that the global 
minimum of f is strictly positive, then we are done. 

Ordinarily, one advantage of a SDP formulation is that global bounds can be found, by consid- 
ering the dual SDP problem, for instance. Thus, the only obstacle remaining is the rank constraint. 
There has been some work done on SDPs with rank constraints, but those are local methods which 
provide low-rank solutions that are not necessarily globally optimal. On the flip side, if we can 
find a low-rank solution at which f attains the value (or if we can simply find a low-rank solu- 
tion to the feasibility problem), then we have found the elusive {5, 3, 3, 3)6 constellation. With this 
in mind, we look at a recently developed and rather controversial method of handling SDPs with 
rank constraints. 

Dattorro has suggested the so-called convex iteration procedure [4] to find low-rank solutions 
to rank-constrained SDPs. Consider a semidefinite program with C as its convex feasible set. Then 
the rank-constrained semidefinite feasibility problem has the form 

find G 

GeS N 

subject to G e C 

(11) 

G>0 

rank(G) < n. 



The feasible set in this case is the intersection of C with a certain subset of the positive semidefinite 
cone boundary, namely the positive matrices of rank n or less. This is clearly a non-convex set, and 
may even be empty. The convex iteration method looks for a feasible solution in this intersection 



7 



by considering the following two coupled SDPs. 



(SDP1) 



minimize tr{Gl¥} 

GsS N 

subject to G e C 



(12) 



G>0, 



and 



minimize tr{GW} 

WeS N 



(SDP2) 



subject to < W < I 



■N 



(13) 



tr{W} =N-n. 

In (fT2~l) . W is an optimal solution to (PT3"1) ; likewise in (PT31 . G is an optimal solution to (PT21 . The 
feasible set in (PT3l) . 



is called the (TV - n)-Fantope. It is the convex hull of the set of all rank-(iV - n) projectors. In fact, 
the extreme points of the (N - n)-Fantope are precisely the rank-(iV - n) projectors. 

One proceeds to solve (PT21) and (PT3l) iteratively, until local convergence of tr{GW} to some 
non-negative value r is established. It is easy to see that the iterations will generate a non- 
increasing sequence of values for tr{GV7}, because if at a certain stage the feasible pair (G', W) 
gives tr {G'VK'} = r', then i J must bound the subsequent minimization problems from above. The 
monotone convergence theorem then guarantees that the sequence of values of tr{GW} converges 
to some real number r > 0. Let us see what we can conclude if we actually have tr {G* W*} = r = 
for some pair (G*, W*). In this case, we would have found a G* e C whose range is orthogonal to 
that of W*. But W* belongs to the (N - n)-Fantope, and has a rank of at least (N - n). Therefore, 
rank(G*) < n, id est we have established feasibility of CCD- 
Returning to our MU constellation existence problem, we see that the convex iteration proce- 
dure allows a systematic way to establish feasibility. We simply construct the rank-constrained 
SDP for the constellation in question, and then use convex iteration to search for a low-rank feasi- 
ble point. An actual numerical implementation readily confirms the existence of the largest known 
constellations such as {5, 3, 3, 2}6 and {5, 5, 5)6, the latter representing three MUBs in C 6 . How- 
ever, convex iteration fails to find (G*, W*) with tr {G* W*} = for the {5, 3, 3, 3} 6 constellation. In 
fact, it appears that the local convergence to r depends on the initialization of the iterative proce- 
dure. The best that has been achieved is r = 0.0022, which is inconclusive. The problem is that 



F 



<N _ 
N-n ~ 



[W eS N :0< W <l N ,iv{W} = N -n), 



(14) 
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it is not known when tr{G*W*} = will actually be achieved via convex iteration for a feasible 
rank-constrained SDR Therefore, our failure to find r = confirms nothing about the existence of 
the {5, 3, 3, 3)6 constellation. Indeed, convex iteration cannot provide an infeasibilty certificate. It 
only adds to the suspicion that a set of four MUBs in C 6 does not exist. 

IV. LASSERRE'S SEMIDEFINITE RELAXATIONS 

Global optimization is extremely difficult for nonconvex problems. Nevertheless, Lasserre has 
developed a remarkable approach to polynomial problems, by defining a sequence of semidefi- 
nite programming relaxations of increasing size which provide ever better approximations to the 
original polynomial problem [^-8]. There is even a publicly available MATLAB implementation 
GloptiPoly 3 [9]. This implementation solves Generalized Problems of Moments (GPM), of which 
the following is a special case: 

minimize L p (x) dp(x) 

^ (15) 
subject to f K hj(x) dp(x) > bj, j = 1, 2, . . . , 

where bj are real numbers and the measure dp in R" is supported on the semialgebraic set K 
defined by the polynomials pu 

K = {x 6 R" : Pi (x) > 0, i = 1, 2, . . .} . (16) 

GloptiPoly 3 carries out the GPM optimization via the moments of the measure dp,, defined as 

y a = [ x a dp(x), oreN", (17) 
Jk 

where a are multi-indices labelling the moments. 

A general constrained polynomial problem has the form 



minimize poW, (18) 
where the set K is defined by the given polynomial constraints as in ([TBI) . Lasserre has shown [6J] 
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that the above polynomial problem can be cast as the moment problem 



minimize f po(x) dp(x) 



subject to f dpi = 1, 



(19) 



i.e., the decision variable dp is a probability measure supported on the set K. A hierarchy of 
semidefinite programming relaxations is then constructed, whose sequence of optimal values con- 
verges to the true global optimal value of (PT9~1) 13,180. 

One can view the equations defining a MU constellation as a system of polynomial equa- 
tions. Following |5|], we consider the simplest constellation that is known not to exist, which 
is {1, 1, 1, 1} 2 . Four real variables suffice to parameterize a candidate set of four vectors: 



(20) 
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A 




+ 1*2, 


' V2 


X-i + L<C4 ; 



Five polynomial equations then serve as constraints defining the MU constellation, 

Pi(x) = x\ + x\ - 1 = 0, 



p 2 (x) = x\ + x\ - 1 = 0, 

p 3 (jc) = (1 + jci) 2 + x\ -2 = 0, 

p 4 (x) = (\+x 3 ) 2 + x 2 4 -2 = 0, 

p 5 (x) = (1 + X1X3 + X2X4) 2 + (X1X4 - X2X3) 2 -2 = 0. 



(21) 



Brierley and Weigert |)5D proposed to take {p\{x)) 2 as the objective function, and solve the opti- 
mization problem 

minimize (p\(x)) 2 

(22) 

subject to pi(x) = 0, i = 2, ... ,5. 

Then, it follows that the global minimum of problem (1221) is zero iff there is a solution to Pi(x) = 
0,i = 1, ... ,5, iff the constellation {1, 1, 1, 1} 2 exists. Now, the optimal objective value attained 
at each relaxation in Lasserre's hierarchy is a global lower bound for the true minimum of the 
original problem. Therefore, a strictly positive optimal objective value at any level of relaxation 
serves as a certificate proving the non-existence of the constellation {1, 1, 1, 1} 2 . In practice, the 
size of the SDP relaxations increases very rapidly as we progress through the hierarchy. In Q5|] , 
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three level of relaxations were calculated for the case of {1, 1, 1, 1} 2 . It was also claimed that the 
first relaxation, which yielded an optimal objective value of 1.4038 x 10~ 8 , provided the desired 
certificate. However, it seems unclear whether this is truly a strictly positive value, or an artifact 
of numerical errors. The second relaxation gives a less ambiguous answer: 0.5359 is certainly 
positive. The authors then proceeded to attempt the same construction for the {5,5,5, 1}6 case, 
which proved to be intractable. In fact, generating the first SDP relaxation was already too difficult, 
let alone solving it. 

The SDP relaxations for constellations like {5, 5, 5, 1}6 and {5, 3, 3, 3}6 are so large because of 
the large number of variables present in the polynomials defining the constellation, as well as 
the relatively high degree of the polynomials, which are quartic. We have tried an alternative 
parameterization of {5, 3, 3, 3}6 by specifying density operators rather than vectors in C 6 . For this, 
we first note that we can choose the density operators for the first group in the constellation to be 
Ef p i = 1, . . . , 5, corresponding the the choice of the computational basis for C 6 . The sixth density 
operator in this group is then automatically determined to be E^ 6 . The remaining 3x3 density 
operators will be parameterized as follows: 

I ZaJA Za,i,2 ZaJ,3 

ZaJA 1 Za,i,(t Za,i,7 

a 1 Za,i,2 Za,i,6 1 Z a ,i,lQ 
Pi ~ ~Z 

Z a ,i,3 ZaJJ Z a ,i,\0 1 
Za,i,4 Z a ,i,% Za,i,\l Z a JA3 
y Za,i,5 Z a ,i,9 Z a ,i,\2 Z a ,i,l4 

The diagonal entries are g because of the MU requirement between the first basis and the remaining 
nine p". Therefore, 15x9 = 135 complex numbers, or 270 real numbers are required. Next, we 
have to ensure that p? is a rank-1 projector. For the rank-1 condition, it is necessary and sufficient 
to check that the columns of p" are linearly dependent, or equivalently, that the determinant of 
every 2nd order minor vanishes. Actually, we can do even better; we only need to impose that 
every neighbouring 2nd order minor vanishes. Since p" is hermitian by parameterization, we do 
not even have to check the minors that involve only the upper triangular entries, since these will 
be duplicated by their complex-conjugate counterparts in the lower triangular sector. Then, the 
projection property follows from the fact that trip"} = 1 is the single non- vanishing eigenvalue 



Za,iA ZaJ,5 

Zaj,8 Z a ,i,9 

Za,i,U Z a ,i,n 

Za,iA3 Za,iA4 

1 Za,iA5 
£<*,/, 15 1 ) 



, a,i = 1,...,3. 



(23) 
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of p". Therefore, the rank-1 projection property is ensured by imposing 15 quadratic equations 
for each p". To enforce the MU constraints between different groups of density operators in the 
constellation, tr{p"p^J = g, a ± j3, i,j = 1,2, 3, an additional 27 quadratic constraints must be 
specified. 

The only remaining constraints are the orthogonality MU constraints 

tr{p^} = d ij , a,i,j= 1,2,3. (24) 

Since we have already ensured that p" is a projector and hence positive semidefinite, we observe 
that these orthogonality MU constraints are simultaneously fulfilled if and only if the single com- 
bined quadratic constraint, 

Po({p-]) = J]K{PiP a j} = 0, (25) 

a 

i<j 

is fulfilled. Hence we can take p to be the objective function to be minimized subject to the 
remaining constraints. This has the nice property that the global minimum of / is zero iff the MU 
constellation {5, 3, 3, 3}(> exists. For this quadratic optimization problem that we have described, 
Gloptipoly 3 is able to generate the first SDP relaxation, which is already huge. It cannot be 
handled on an ordinary desktop PC. Although this is a slight improvement, it is likely that this first 
relaxation will yield inconclusive results, if the smaller-sized problems are to be a guide. 

It is worth noting that reducing the problem from a quartic one to a quadratic one comes at the 
cost of introducing more variables and constraint equations. On a positive note, QCQP is itself an 
active field of research, so there is some hope that the large QCQP problem that we have described 
above can actually be handled. 

V. ALGEBRAIC GEOMETRY AND GROBNER BASES 

It might very well be that carrying out a polynomial optimization to prove the non-existence 
of a certain MU constellation is an overkill. Perhaps it is really unnecessary to optimize or solve 
the polynomial equations; we might be satisfied with just knowing how the solutions "look like", 
or how many of them there are. In essence, what we have is a set of N multivariate polynomial 
equations {pi(x) = 0} i= i.. jv defining a MU constellation, and what we are interested in are the 
solutions, if any, to these equations. This leads us to the field of algebraic geometry, which is 
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replete both with elegant theorems and frustrating open problems. 

The polynomials p, belong the the ring of polynomials over R in variables x\,..., x n , which is 
denoted by R[xi, . . . , x n ]. Let S be the set [pi}i=i,...,jv- Then the central object of interest is the set 
of common zeros of S in R", i.e. the (real) variety Vr(S) = {x e R" : p(x) = Vp e S). In this 
definition, the MU constellation defined by S exists if and only if Vr(S) + 0. Now, real algebraic 
geometry is notoriously difficult, so we look for complex zeros instead. That is, we consider S 
as a subset of C[xi, . . . , x n ] and study the complex variety V C (S) = {x e C" : p(x) = Vp e S }. 
Certainly, if V C (S) = 0, then V K (S) = 0. 

Different sets of polynomials can give rise to the same algebraic variety. For example, if x is a 
zero of polynomials / and g, then it is also a zero of f+g and fg. Therefore, Vc(S) is equivalently 
given by V c ((pi,. ■ -,Pn)), where (pi,.. .,p N ) is the ideal generated by S ={pi,.. .,p N ], 

(p u ...,p N ) = jpe C[xi,...,x n ] : p = ^ r iPh n e C[x l5 . . . ,x n ], i = l,...,ivj. (26) 

The ideal (p lt . . . , p N ) is not uniquely generated by S . In fact, every ideal / in C[xi, . . . , x n ] is 
finitely generated, and if / is non-zero, there even exists a unique, distinguished generating set, the 
reduced Grobner basis (w.r.t. a monomial order, to be defined later), which generates the ideal. 
Grobner bases can be viewed as a generalization of Gaussian elimination for linear systems or 



the Euclidean algorithm for computing greatest common divisors (see, for example, HOD). The 
notion of a Grobner basis only makes sense after one has defined some monomial order, i.e., a 
total order on the set of all monic polynomials in a polynomial ring which respects multiplication 
(u < v => uw < vw), and is a well-ordering (every non-empty set of monomials has a minimal 

element). For a single variable, the only monomial ordering is 1 < x < x 2 < x 3 For the 

general multivariate case, one example of a monomial order is the lexicographic order (lex). This 
firstly requires that x x > x 2 > ■ ■ ■ > x n . For higher degree monomials, the exponents of X\ are 
compared; in the event of a tie, the exponents of x 2 are compared, and so on. Thus, for instance, 
XjX4 > x\xl > x\x 2 x\. Other examples of monomial orders are the graded lexicographic order 
(grlex) and the graded reverse lexicographic order (grevlex). In grlex, the total degree of the 
monomials are first compared, and ties are broken by applying lex. In grevlex, total degree is 
first compared; ties are broken by comparing exponents of x„, with smaller exponents regarded as 
larger in the ordering, followed, if necessary, by comparing exponents of x n ^,x n - 2 , etc. Given a 
monomial order and a polynomial p 6 C[xi, . . . , x„], we denote the largest monomial in p by lp(p), 
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and its corresponding coefficient by lc(p). With respect to a given monomial order, a Grobner basis 
is defined as follows: 

Definition 1. A set of non-zero polynomials G = {gi, ■ ■ ■ ,g t ) contained in an ideal I, is called a 
Grobner basis for I iff for all non-zero f £ I, there exists i £ {1, . . . ,t) such that lp(g,) divides 

W)- 

A Grobner basis for an ideal is not unique, but the reduced Grobner basis is. This is defined as 
follows. 

Definition 2. A Grobner basis G = {gi,. . .,g t } is called a reduced Grobner basis if for all i, 
lc(g,-) = 1, and no non-zero monomial in g,- is in the ideal generated by the leading terms of the 
elements in G - {gi}. 

The importance of Grobner bases lies in the following theorem: 

Theorem 1. Let I be an ideal ofR[x\, . . . , x n ], and G be its reduced Grobner basis with respect to 
a monomial order. Then, VcCD = if and only if I £ G, or equivalently, if and only ifG = {1}. 

Therefore, the system of real polynomial equations p\ = 0, P2 = 0, . . . , Pn = has no solutions 
in C' ! iff the reduced Grobner basis for the ideal (pi, . . . ,pn) is equal to {1}. In [5], the authors 
confirmed that the constellation {1, 1, 1, 1)2 does not exist by computing that the polynomial equa- 
tions defining the constellation have G = {1} as the reduced Grobner basis. Thus, we see that we 
can find out something about the solution set of a system of polynomials (its non-existence for 
instance), without actually finding the solutions directly. This is the advantage of algebraic geo- 
metric methods such as Grobner basis computations. Of course, we should have a way of finding 



Grobner 



5ases in the first place. Fortunately, there is a very general algorithm, Buchberger's Al- 



gorithm [11], that can compute the Grobner basis for the ideal I = {p\, . . .,Pn) in a finite number 
of steps. Unfortunately, algorithms for computing Grobner bases are not as straightforward as, 
say, Gaussian elimination. In general, Grobner bases can be very large, and the computational 
cost of finding one depends very much on the monomial order, the order of the polynomials, and 
the choice of the so-called S -polynomials which appear in the intermediate steps of Buchberger's 
Algorithm. For instance, it was reported in Q3D that the computation of a Grobner basis for the MU 
constellations {5, 3, 3, 3}6 and {5, 5, 4, 1} using the package FGb failed because of memory issues 
(despite having 16GB of memory). Grobner basis computation is still actively researched, and 



other algorithms such as Faugere's F4 Ill2h and F5 II 1311 algorithms are available 
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VI. HILBERT'S NULLSTELLENSATZ, NULLA, AND PARRILO'S SUM-OF-SQUARES 

Computation of a Grobner basis may provide an infeasibility certificate for a system of poly- 
nomials. However, even that might be more than what we require. There is a beautiful theorem by 
David Hilbert, which in one particular form states the following: 

Theorem 2 (Hilbert's Nullstellensatz). A system of polynomial equations, p\{x) = 0, . . . ,p s {x) = 
0, has no solutions over an algebraically closed field K if and only if there exist polynomials 
ri,,..,r s e K[x\ ,...,x n ] such that 1 = £;=i npu 

Therefore, an identity of the form 1 = 2/=i r >P> provides an infeasibilty certificate for the 
system p\(x) = 0, . . .,p s (x) = 0. The maximum degree of the polynomials ript will be called the 
degree of the Nullstellensatz certificate. We remark that the Nullstellensatz is intimately linked 
to Theorem [T]on Grobner bases. We also know that a Nullstellensatz certificate must exist for an 
infeasible system of polynomials. Therefore, we look to find such a Nullstellensatz refutation for 
the existence of certain MU constellations. There are in fact systematic ways to search for such 



15] 



refutations, for instance, NulLA (Nullstellensatz Linear Algebra) [14, 

The basic idea behind NulLA is quite simple. We fix a tentative degree d for the Nullstellensatz 
certificate that we wish to find. We then expand the assumed polynomial identity 1 = Yfi=\ r tPi mt0 
a linear combination of monomials with degrees less than or equal to d. The coefficients of these 
monomials will be linear expressions in the coefficients defining the unknown polynomials r,-. 
Note that two polynomials over a field are equal if and only if the coefficients of every monomial 
are equal. Therefore, the identity 1 = 2;=i r iPi corresponds to a linear system of equations in 
the coefficients of r ; . Solving this linear system then results in two possible outcomes. If the 
system is consistent, then any solution produces the desired Nullstellensatz certificate of degree 
d. Otherwise, no certificate of degree less than or equal to d exists, and we start afresh with a 
tentative degree d + 1. Repeat the process if necessary, until a certificate is found. 

The linear systems that appear in NulLA increase very rapidly in size as d increases, and can 

be huge for reasonably- sized problems, even for d as small as six. We can see this easily. There 
In + d\ 

are I I monomials in n variables of degree d or less. Writing dj = d - deg(p,), there are 

2jU | j.^ j unknowns in the linear system to be solved for a Nullstellensatz certificate of degree 
d. Upper bounds on the degree of a Nullstellensatz certificate are known to be doubly-exponential 



in the number of input polynomials and their degree. However, as pointed out in [11411 . fairly 
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low-degree certificates exist for many examples. 

To see how NulLA works out in practice, we look once again at the constellation {1, 1, 1, 1} 2 . 
The equations defining it are given in (|2"TT) . After carrying out the NulLA algorithm, we find a 
Nullstellensatz certificate of degree 6, given by 



n(x) = 




— 2jc^ + X\x\ + X2X3X4. 


2 3 

— X2X 1) X4 — X%X^ 


r 2 (x) = 




- X3 + 2x^ — X2X4J , 




r 3 0) = 




2 \ 
Xi-tj — X2X3X4 1 , 




r 4 (x) = 








rs(x) = 


2X2X4- 







(27) 



It is straightforward to verify tha t r,/?* = 1. We should point out that during the Grobner 
basis computation carried out in I5D, a slightly different degree 6 Nullstellensatz certificate was 



produced as a by-product. 

Suppose we are looking for a certificate of degree d. Then the linear system we need to solve 

(n + d\ 

looks like My = b. Here, the matrix M has I I rows, one per monomial x a of degree d or 

less (or is a multi-index); M also has one column per polynomial of the form x 6 f u where x 5 is a 

In + d\ 

monomial of degree less than or equal to d - deg(/i). The vector b has entries which are 



d 

zero everywhere except for the entry corresponding to the constant monomial x°, which is 1 . From 
this description, one sees that the size the matrix M grows very quickly with the certificate degree 
d. 



Still, there are a number of ways to optimize NulLA II 1411 . Firstly, we note that the system 
of polynomials p y , . . . , p s defining a MU constellation has only real coefficients. Then it can be 
shown that it suffices to search for real Nullstellensatz certificates. In other words, a Nullstellensatz 
certificate 1 = £ • r,/>, where r, e C[x\ . . . , x n ] exists iff there exists a real Nullstellensatz certificate 
1 = Yi S i r [Pi where r\ e R[jci, . . . ,x n ]. Linear algebra can then be carried out over the reals. 
Secondly, the size of the linear system My = b that has to be solved at each stage of NulLA can be 
significantly reduced if there are certain symmetries in the system of equations p u . . . , p s = 0. For 
instance, suppose that the set S = {p, } != i..., s is invariant under the action of a group of permutations 
G of the variables x\,...,x n . We denote the image of p { under g e G by g(fi). G also induces 
an action on the set of monomials of degree t. The orbit of a monimial x a under G is denoted by 
0{x"), while the orbit of x 6 fi is denoted by Oix 6 fi). In view of the symmetries captured by G, we 
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introduce the matrix equation My = b, where the rows the matrices are indexed by the orbits 0{x") 
and the columns indexed by the orbits 0(x s fi). The entries of M are defined by 



M, 



(28) 



Note that this definition is independent of the choice of x a in the orbit Oix"). The vector b has 
zeroes everywhere except for the entry corresponding to 0(1). With these definitions, it can be 
shown that if a solution to My = b exists, then a solution to My = b, exists. 

The above idea may be applicable when dealing with, for instance, the {5, 5, 5, 1 constellation. 
In this case, we may choose the first group of vectors to be the computational basis vectors, and the 
last singleton to be 4=(1, 1, 1, 1, 1, The remaining two groups of five vectors can be written as 





f 


1 ' 






Xi 






1 


Xi 


7,2 + iyua 




V6 


Xi 


7,3 + iji.7,3 


> 




Xi 


7,4 + ty,JA 






\Xi 


7,5 + %7;5, 





(29) 



/=1,2;;=1,...,5 



with the index i labelling the group of vectors, and j labelling the five vectors in each group. In 
this parameterization, there are 100 real variables, Xij t k,yij,k, with k e {1,2,3,4,5} labelling the 
component number. The MU equations are: 



2 -1=0 



(l+ZLi *«,/,*) +(Xa ,.v/.a) "6 = 

(l + ZLi (xij, k Xi.j, k + yij tk yi,j tk )f + (xij,kyi>,f,k ~ xvj&UjjSf = 



Vi, j, k, 

V/,7, 
6 Vi*i',Vj,f 
Mi = i',\fj±f 



(30) 



It follows that the above set of equations are invariant under simultaneuous permutations 

{xi,j,k\y?,r,k>} -> {^(o^w^wJ^v^y^of)) in each of the three indices ' wnere g 1 ,^,^ are > 

respectively, permutations from the permutation groups on 2, 5, and 5 objects. In view of this, the 
linear system My = b should be significantly simpler than the original system My = b. 



A number of other possible ways to improve the efficiency of NulLA are described in II14 . 
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150. These include appending extra polynomial equations from the radical ideal of (pi, . . . , p s ) to 
reduce the Nullstellensatz degree, branching the polynomial system into smaller subsystems with 
the aim of finding lower-degree infeasibility certificates for the smaller subsystems, and using 
alternative versions of the Nullstellensatz itself. A recent modification to the NulLA algorithm, 



called FPNulLA (fixed-point NulLA), is also proposed in [16]. 

So far, we have been trying to show that there are no common complex zeros of the polynomial 
system S defining a certain MU constellation. It could well be the case that V R (S) = but V C (S) # 
0. In other words, a MU constellation may be non-existent, but no Nullstellensatz certificate can 
prove it. Unfortunately, the real Nullstellensatz is somewhat more complicated than Hilbert's 
original one. One version says the following: 

Theorem 3 (Real Nullstellensatz). Let p\,...,p s £ R[x\ , . . . , x n ], and let Z c R[jci , . . . , x n ] denote 
the cone of polynomials representable as a sum-of -squares. Then Vr((pi, . . . , p s }) = if and only 
if -I eS + (pi,...,p s ). 



The real Nullstellensatz is closely related to Stengle's Positivstellensatz Ill7h . which deals with 
the semialgebraic set 

S S A = {x £ R" : Pi(x) = and qj(x) > Mi = 1, . . . , s; j = 1, . . . ,t] , (31) 

where p t , qj £ R[xi, . . . , x n ]. We define M({qj}) to be the set of all finite products of qj, including 
the empty product, 1. The cone associated with [qj] is defined as 

cone(q u ...,q t ) = \go + ^gkb k ■ go, ■ ■ ■ ,g r 6 Z,b u . . . ,b r £ M({qj}» (32) 



The Positivstellensatz reads 

Theorem 4 (Positivstellensatz). The semialgebraic set S sa = if cind only if -I £ 
cone(qu...,q t ) + (pu.--,p s ). 



In 111811 . Parrilo described how the search for a Positivstellensatz infeasibility certificate can be 
rephrased as a hierarchy of semidefinite programs. This essentially involves assuming a certain 
maximal degree for the infeasibility certificate, and checking the feasbility of the corresponding 
SDP If the SDP for a certificate of degree d is not feasible, then one proceeds to check the feasi- 
bility of the SDP for a d + 1 degree certificate. 



18 



Direct numerical searches have not been able to find MU constellations such as {5,3,3,3)6 
and {5,5,4, 1} 6 . If such constellations were to exist, it is possible that the corresponding ideal 
S = (pi, . . . ,p s ) is zero-dimensional, i.e., the set of common zeros, Vr(5), is finite. In this case, 
Lasserre et. al. J19I1 have described a numerical algorithm based on semidefinite programming to 
compute the points on this finite variety. Finally, we remark that infeasibility certificates could be 
easier to find for systems that are "more infeasible". From this point of view, one might want to 
consider directly the case of seven MUBs in C 6 instead of MU constellations. 



VII. CONCLUSION 

A few approaches to the existence problem of four MUBs in C 6 have been described. Current 
numerical evidence suggests that four MUBs do not exist, but a rigorous infeasibility certificate 
is lacking. One approach to searching for such an infeasibility certificate involves setting up an 
optimization problem in which the objective function attains the maximum/minimum value in 
the codomain precisely when a set of four MUBs exists. Then, one can obtain an infeasibility 
certificate by proving certain global bounds on this objective function. The main obstacle in this 
approach is the absence of convexity, which appears to be crucial; relaxation of the non-convex 
constraints seems to have the tendency to give trivial global bounds. 

A different approach views the problem from the point of view of algebraic geometry. A MU 
constellation is defined using a set of multivariate polynomial equations, and its non-existence 
corresponds to the infeasibility of this set of equations. A number of computable infeasibilty 
certificates are available from the theory of algebraic geometry, but these usually prove the non- 
existence of complex zeros. Furthermore, the number of variables and constraints required to 
define a MU constellation is quite large, and poses problems for realistic computation of infeasi- 
bilty certificates. One recent idea is to use linear algebra to compute infeasibilty certificates, and it 
is hoped that with the appropriate refining techniques, the linear systems involved can be solved. 
Real algebraic geometry is less well-understood, but recent work has linked it to semidefinite pro- 
gramming, which may allow infeasibilty certificates to be computed with a realistic amount of 
resources. 

Special thanks to Philippe Raynal for numerous valuable discussions, and B.-G. Englert for his 
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